
function [Y,f]=mydft(y,Fs,b)
%[Y,f]=mydft(y,Fs)
%y: data vector
%Fs: sampling frequency
%Y:complex fft value
%f: frequency data vector
%b: 0 single-sided spectrum
%   1: double-sided spectrum

L=length(y);
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
if b==0
    f = Fs/2*linspace(0,1,NFFT/2);
else
    f = -NFFT/2:NFFT/2-1
    f=f*Fs/NFFT;
    %Fs/2*linspace(-1,1,NFFT);
end
if b==0
    % Plot single-sided amplitude spectrum.
    plot(f,2*abs(Y(1:NFFT/2))) 
    title('Single-Sided Amplitude Spectrum of y(t)')
else
    % Plot double-sided amplitude spectrum.
    %plot(f(1:end),[abs(Y((NFFT/2)+1:end))  abs(Y(1:NFFT/2))]) 
    plot(f,abs(fftshift(Y)));
    title('Double-Sided Amplitude Spectrum of y(t)')
end
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
end